logo of company

From Risk to Resilience: Ecological Exploration of Psychopathology and Flourishing in Neurodivergent Children Facing Adversity


This study extends analyses conducted by Bethell et al. (2023), who explored the role of family resilience and connection in promoting flourishing among U.S. children facing adversity. Our study applies a similar framework but focuses on neurodivergent children, examining how adverse childhood experiences (ACEs) influence mental health outcomes and how flourishing behaviors—such as curiosity, emotional control, and task persistence—may buffer these effects. This study was conducted by the Slopen Laboratory at the Harvard T.H. Chan School of Public Health.

Author: Adrián Medina

Date: December 5, 2024

Project Overview


This study investigates the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.

Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 54,528, MAge = 12.1). This cross-sectional analysis assesses the exposure to ACEs and its impact on three key mental health outcomes: anxiety, depression, and behavioral problems. Interaction terms in logistic regression models are used to assess whether child flourishing behaviors buffer the effects of ACEs on these outcomes.

Data Collection:
Key variables collected include:

  • Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
  • Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
  • Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).

Analytic Approach:
Logistic regression models were employed to evaluate the association between ACEs and mental health outcomes, including interaction terms to investigate the moderating role of child flourishing. The models incorporated a comprehensive set of covariates, including child age, sex, race/ethnicity, socioeconomic status, caregiver mental health, neighborhood conditions, and healthcare access metrics.

Goals:
The study aims to:

  • Examine the dose-response relationship between ACE exposure and mental health outcomes among neurodivergent children.
  • Assess the protective role of child flourishing behaviors, with the goal of identifying resilience-promoting factors that can inform interventions for this vulnerable population.
Tip

For detailed information on the NSCH-derived variables, refer to the NSCH Data Dictionary.

For a primer on Neurodiversity: What is Neurodiversity?

Data Frame Initialization


Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.

Expand Code
# Set CRAN repository for consistent package installation
options(repos = c(CRAN = "http://cran.rstudio.com/"))

# Load or install the pacman package for efficient package management
if (!require(pacman)) install.packages("pacman")

# Use p_load to load or install the required packages
p_load(
  dplyr,           # Data manipulation
  tidyr,           # Data reshaping (e.g., pivot_longer)
  readr,           # Reading data (CSV files)
  knitr,           # Creating tables for reports
  kableExtra,      # Enhancing `knitr::kable()` tables
  ggplot2,         # Data visualization
  sjlabelled,      # Labelled data utility functions
  sjmisc,          # Additional data manipulation functions
  survey,          # Complex survey analysis
  lmtest,          # Testing models (Breusch-Pagan test) and robust standard errors
  sandwich,        # Calculating robust standard errors
  gtsummary,       # Creating regression summary tables
  ggeffects,       # Creating prediction plots
  gt,              # Adjusting frequency/contingency tables
  psych,           # Factor analysis
  crosstable,      # Cross-Tabulations
  ggpubr,          # Summary statistics + visualizations
  tidyverse,       # Data visualization supplement
  ggdist           # Distribution visualizations
  )

# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WD
base_path <- "~/GitHub/Adversity-Neurodiversity/data_files"
setwd(base_path)

# Load primary data file
neurodivergent_data <- read_csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2016to2023.csv")
1
See details under the Data Extraction & Filtering (Archived Code) callout below!
Data Extraction & Filtering (Archived Code)

The file being used for these analyses is a subset of a “master” file, Data2_2016to2022.csv, containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is >500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.

Additionally, data from the NSCH 2023 was also later merged into the subset file being used for this analyses. The code for the extraction/wrangling process is shown below.

Expand Code
# Define the path to the dataset
data_path <- "~/Downloads/Data2_2016to2022.csv"  
Data2_2016to2022 <- read.csv(data_path)

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status
neurodivergent_data <- Data2_2016to2022 %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis,
         k2q35a, k2q35c, k2q38a, k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, 
         k2q37a, k2q37c, downsyn, downsyn_desc, k2q31a, k2q31c, k2q61a, cerpals_desc, 
         k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, 
         k7q85_r, k7q84_r, talkabout, wktosolve, strengths, hopeful, k8q21, k8q30, fwc, 
         k2q42a, k2q42c, a1_menthealth, a2_menthealth, k2q01, k4q27, k4q28x01, k4q28x04, 
         appointment, available, issuecost, notelig, notopen, transportcc, family, family_r,
         k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k10q40_r, k10q41_r,
         k10q30, k10q31, k9q96, sc_cshcn, k3q22, k3q20, instype, treatneed, k4q22_r,
         menbevcov) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1 | k2q42a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data, "neurodivergent_data_2016to2022.csv")
Expand Code
# Define the path to the dataset
data_path <- "~/Downloads/nsch_2023e_topical_SAS/nsch_2023e_topical.csv"  
nsch_data <- read.csv(data_path)

# Create a new variable 'fpl_mean' representing the average of specific columns
nsch_data <- nsch_data %>%
  mutate(fpl_mean = rowMeans(select(., fpl_i1, fpl_i2, fpl_i3, fpl_i4, fpl_i5, fpl_i6), na.rm = TRUE))

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status (manually added in blank columns fpl and downsyn_desc)
# `menbevcov` was not collected for 2023
neurodivergent_data_23 <- nsch_data %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, 
         k2q35a, k2q35c, k2q38a, k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, 
         k2q37a, k2q37c, downsyn, downsyn_desc, k2q31a, k2q31c, k2q61a, cerpals_desc, 
         k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, 
         k7q85_r, k7q84_r, talkabout, wktosolve, strengths, hopeful, k8q21, k8q30, fwc, 
         k2q42a, k2q42c, a1_menthealth, a2_menthealth, k2q01, k4q27, k4q28x01, k4q28x04, 
         appointment, available, issuecost, notelig, notopen, transportcc, family_r,
         k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k10q40_r, k10q41_r,
         k10q30, k10q31, k9q96, sc_cshcn, k3q22, k3q20, instype, treatneed, k4q22_r) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1 | k2q42a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data_23, "neurodivergent_data_2023.csv")

## Manually removed index column in both 2016-2022 and 2023 files 
## Manually add empty `family` & `menbevcov` columns to 2023 file
# Read in update 2023 data file
neurodivergent_data_23 <- read.csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2023.csv")

# Read in prior years of neurodivergent data
neurodivergent_data_16to22 <- read.csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2016to2022.csv")

# Merge two data frames with identical columns
neurodivergent_data_16to23 <- rbind(neurodivergent_data_16to22, neurodivergent_data_23)

# Write the merged data frame into a CSV file
write.csv(neurodivergent_data_16to23, "neurodivergent_data_2016to2023.csv", row.names = FALSE)

Analytic Data Preparation


Child Demographic Variables

Recode survey year into a factor, recode child sex as a categorical variable, recode race and ethnicity into descriptive categories, and create a combined race-ethnicity variable.

Expand Code
# Recode survey year into a factor variable
neurodivergent_data <- neurodivergent_data %>%
  mutate(survey_year = factor(year))

# Recode sex and race/ethnicity categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode sex into categorical factors
    sc_sex_cat = factor(sc_sex, levels = c(1, 2), labels = c("Male", "Female")),
    
    # Recode race categories into descriptive labels
    race = case_when(
      sc_race_r == 1 ~ "White",
      sc_race_r == 2 ~ "Black",
      sc_race_r == 3 ~ "American Indian or Alaska Native",
      sc_race_r %in% 4:5 ~ "Asian, Native Hawaiian, or Pacific Islander",
      sc_race_r %in% 6:7 ~ "Other or Mixed Race",
      TRUE ~ NA_character_
    ),
    
    # Recode ethnicity
    ethnicity = case_when(
      sc_hispanic_r == 1 ~ "Hispanic/Latino",
      sc_hispanic_r == 2 ~ "Non-Hispanic/Latino",
      TRUE ~ NA_character_
    ),
    
    # Create combined race and ethnicity category
    sc_race_cat = case_when(
      race == "White" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic White",
      race == "Black" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Black or African American",
      race == "American Indian or Alaska Native" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic American Indian or Alaska Native",
      race == "Asian, Native Hawaiian, or Pacific Islander" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander",
      race == "Other or Mixed Race" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Other or Mixed Race",
      ethnicity == "Hispanic/Latino" ~ "Hispanic/Latino of Any Race",
      TRUE ~ NA_character_
    ),
    sc_race_cat = factor(sc_race_cat, levels = c(
      "Non-Hispanic White", 
      "Non-Hispanic Black or African American", 
      "Non-Hispanic American Indian or Alaska Native", 
      "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", 
      "Non-Hispanic Other or Mixed Race", 
      "Hispanic/Latino of Any Race"
    ))
  )

Family & Caregiver Context Variables

Recode caregiver education, family income, and caregiver composition variables to categorical values. Create a composite score for caregiver mental health, then categorize it into meaningful health levels.

Expand Code
# Recode caregiver education and family income categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode caregiver education into categories
    caregiver_education = case_when(
      higrade_tvis == 1 ~ "Less Than High School",
      higrade_tvis == 2 ~ "High School/Vocational",
      higrade_tvis == 3 ~ "Some College/Associate Degree",
      higrade_tvis == 4 ~ "College Degree/Higher",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    caregiver_education = factor(caregiver_education, levels = c("Less Than High School", "High School/Vocational", "Some College/Associate Degree", "College Degree/Higher")),

    # Recode family income into federal poverty level categories
    fpl_category = case_when(
      fpl %in% c(50:99) | fpl_mean < 100 ~ "<100%",
      fpl %in% c(100:199) | fpl_mean < 200 ~ "100%-199%",
      fpl %in% c(200:299) | fpl_mean < 300 ~ "200%-299%",
      fpl %in% c(300:399) | fpl_mean < 400 ~ "300%-399%",
      fpl %in% c(400:999) | fpl_mean < 999 ~ "≥400%",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    fpl_category = factor(fpl_category, levels = c("<100%", "100%-199%", "200%-299%", "300%-399%", "≥400%"))
  )

# Simplify caregiver composition based on family structure
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    caregiver_composition = case_when(
      family %in% c(1, 2, 3, 4) | family_r %in% c(1, 2, 3, 4) ~ "Dual-caregiver household",
      family %in% c(5, 6, 7, 8, 9) | family_r %in% c(5, 6, 7, 8) ~ "Single-caregiver household",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    caregiver_composition = factor(caregiver_composition, levels = c("Dual-caregiver household", "Single-caregiver household"))
  )

# Create a composite variable for caregiver mental health
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    caregiver_mental_health_score = rowMeans(select(., a1_menthealth, a2_menthealth), na.rm = TRUE),
    caregiver_mental_health = cut(
      caregiver_mental_health_score,
      breaks = c(1, 1.5, 2.5, 3.5, 4.5, 5),
      labels = c("Excellent", "Very Good", "Good", "Fair", "Poor"),
      include.lowest = TRUE
    )
  )
1
Scores of caregiver mental health are averaged to account for single respondents (i.e., single-caregiver households).

Neighborhood Context Variables

Recode neighborhood binary items to 0/1 format and standardize ordinal neighborhood variables by subtracting their mean and dividing by the standard deviation.

Expand Code
# Recode neighborhood condition-related binary items to 0/1
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    across(c(k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k9q96), ~ ifelse(. == 1, 1, 0))
  )

# Standardize neighborhood condition-related ordinal variables
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    across(c(k10q40_r, k10q41_r, k10q30, k10q31), ~ scale(., center = TRUE, scale = TRUE))
  )
1
k10q11 = Neighborhood - Sidewalks or Walking Paths; k10q12 = Neighborhood - Park or Playground; k10q13 = Neighborhood - Recreation Center; k10q14 = Neighborhood - Library or Bookmobile; k10q20 = Neighborhood - Litter or Garbage; k10q22 = Neighborhood - Poorly Kept or Rundown Housing; k10q23 = Neighborhood - Vandalism; k9q96 = Other Adult Child Can Rely On For Advice
2
k10q40_r = Child is Safe In Neighborhood; k10q41_r = Child Is Safe at School; k10q30 = People In Neighborhood Help Each Other Out; k10q31 = Watch Out for Other’s Children

Factor Analysis - Neighborhood Context Scores

Conduct factor analysis on neighborhood context variables to derive factor scores for physical resources, environmental threats, and collective efficacy, and add these scores to the primary data frame. Remove intermediate data frames from the global environment.

The factor structure was adapted from factor conceptualization detailed in Fan & Chen (2012), “Family functioning as a mediator between neighborhood conditions and children’s health: Evidence from a national survey in the United States.” Principal Axis Factoring (PAF) was used as the factoring method, and factor scores were calculated using the regression method, thus consistent with the approach in the referenced study.

Expand Code
# Select the items for each neighborhood factor
physical_resources_items <- neurodivergent_data %>%
  select(k10q11, k10q12, k10q13, k10q14)

environmental_threats_items <- neurodivergent_data %>%
  select(k10q20, k10q22, k10q23, k10q40_r, k10q41_r)

collective_efficacy_items <- neurodivergent_data %>%
  select(k10q40_r, k10q41_r, k10q30, k10q31, k9q96)

# Conduct factor analysis for each factor using Principal Axis Factoring (PAF)
fa_physical_resources <- fa(physical_resources_items, nfactors = 1, fm = "pa")
fa_environmental_threats <- fa(environmental_threats_items, nfactors = 1, fm = "pa")
fa_collective_efficacy <- fa(collective_efficacy_items, nfactors = 1, fm = "pa")

# Calculate factor scores for physical resources using regression method
physical_resources_scores <- factor.scores(physical_resources_items, fa_physical_resources, method = "regression")$scores[, 1]

# Calculate factor scores for environmental threats using regression method
environmental_threats_scores <- factor.scores(environmental_threats_items, fa_environmental_threats, method = "regression")$scores[, 1]

# Calculate factor scores for collective efficacy using regression method
collective_efficacy_scores <- factor.scores(collective_efficacy_items, fa_collective_efficacy, method = "regression")$scores[, 1]

# Add factor scores back to the dataset
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    physical_resources = physical_resources_scores,
    environmental_threats = environmental_threats_scores,
    collective_efficacy = collective_efficacy_scores
  )

# Remove intermediate data frames from the global environment
rm(physical_resources_items, environmental_threats_items, collective_efficacy_items, 
   fa_physical_resources, fa_environmental_threats, fa_collective_efficacy, 
   physical_resources_scores, environmental_threats_scores, collective_efficacy_scores)

Child Health & Development Variables

Standardize Down Syndrome & Tourette Syndrome severity variables based on survey year, calculate total severity scores across diagnoses, categorize average severity scores into severity levels, recode child general health, and convert the special healthcare needs indicator into binary format with labeled categories.

Expand Code
# Standardize Down Syndrome and Tourette Syndrome severity variables based on survey year
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    downsyn_desc_recode = case_when(
      survey_year %in% c(2016, 2017) & downsyn_desc == 1 ~ 1,
      survey_year %in% c(2016, 2017) & downsyn_desc == 4 ~ 3,
      survey_year >= 2019 ~ downsyn_desc,
      TRUE ~ NA_real_  # Handle missing or undefined years as NA
    ),
    k2q38c_recode = case_when(
      k2q38c == 1 ~ 1,
      k2q38c == 4 ~ 3
    )
  )

# List of severity variables, including recoded Down & Tourette Syndromes' description
severity_variables <- c("k2q31c", "k2q35c", "k2q38c_recode", "k2q30c", "k2q36c", "k2q60c", 
                        "k2q37c", "downsyn_desc_recode", "cerpals_desc", "k2q42c")

# Calculate total and average severity scores
neurodivergent_data <- neurodivergent_data %>%
  rowwise() %>%
  mutate(
    total_severity_score = ifelse(all(is.na(c_across(all_of(severity_variables)))), NA_real_, 
                                  sum(c_across(all_of(severity_variables)), na.rm = TRUE))
  ) %>%
  ungroup()

# Recode child general health variable
neurodivergent_data <- neurodivergent_data %>%
  mutate(child_gen_health = factor(
    case_when(
      k2q01 == 1 ~ "Excellent",
      k2q01 == 2 ~ "Very Good",
      k2q01 == 3 ~ "Good",
      k2q01 == 4 ~ "Fair",
      k2q01 == 5 ~ "Poor"
    ),
    levels = c("Excellent", "Very Good", "Good", "Fair", "Poor"))
    )

# Convert Special Healthcare Needs (SHCN) indicator to binary and label accordingly
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    sc_cshcn = ifelse(sc_cshcn == 1, 1, 0),
    special_healthcare_needs = factor(
      case_when(
        sc_cshcn == 0 ~ "Non-SHCN",
        sc_cshcn == 1 ~ "SHCN"
      ),
      levels = c("Non-SHCN", "SHCN")
    )
  )
1
k2q31c = ADD/ADHD Severity Description; k2q35c = Autism ASD Severity Description; k2q38c_recode = Tourette Syndrome Severity Description; k2q30c = Learning Disability Severity Description; k2q36c = Developmental Delay Severity Description; k2q60c = Intellectual Disability Severity Description; k2q37c = Speech Disorder Severity Description; downsyn_desc_recode = Down Syndrome Severity Description; cerpals_desc = Cerebral Palsy Severity Description; k2q42c = Epilepsy Severity Description

Psychosocial Risk & Protective Factors

Recode ACEs, psychopathology, and psychosocial factors. Compute the Childhood Flourishing Index and the Family Resilience and Connection Index, based on Bethell et al. (2019), “Family Resilience And Connection Promote Flourishing Among US Children, Even Amid Adversity.” Create categorical versions of these indices for further analysis.

  • The Child Flourishing Index (CFI) in was calculated using three items from the National Survey of Children’s Health (NSCH). Each item asked parents how well each of the following statements described their child: (1) “shows interest and curiosity in learning new things,” (2) “works to finish tasks he or she starts,” and (3) “stays calm and in control when faced with a challenge.” Each item was assigned one point for a response of “definitely true.” Children who scored 3 were classified as flourishing.

  • The Family Resilience and Connection Index (FRCI) was created using six items from the NSCH. The index included four items on family resilience, which asked how often families: (1) “talk together about what to do,” (2) “work together to solve our problems,” (3) “know we have strengths to draw on,” and (4) “stay hopeful even in difficult times.” Parents could earn one point for each of these items if they responded “all of the time.” Two additional items assessed family connection: (1) the ability of the parent to “share ideas or talk about things that really matter” with their child and (2) how well parents were “handling the day-to-day demands of raising children.” One point was assigned for each response of “very well.” The FRCI score ranged from 0 to 6.

Expand Code
# Recode Adverse Childhood Experiences (ACE) variables & psychopathology outcomes
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode 'ace1' to a dichotomous variable
    ace1_recode = ifelse(ace1 %in% c(2, 3, 4), 1, ifelse(ace1 == 1, 0, NA))
  ) %>%
  mutate(
    # Total ACE count excluding missing values
    ace_total = rowSums(select(., ace1_recode, ace3:ace12) == 1, na.rm = TRUE),
    
    # Categorize total ACEs for analysis
    ace_counts = cut(ace_total, breaks = c(-1, 0, 1, 3, Inf), labels = c('0 ACEs', '1 ACE', '2-3 ACEs', '4+ ACEs'), right = TRUE),
    ace_counts = factor(ace_counts, levels = c("0 ACEs", "1 ACE", "2-3 ACEs", "4+ ACEs")),
    
    # Recode psychopathology outcomes
    Anxiety = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    Depression = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    Behavioral_Problems = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes"))
  )

# Calculate and recode Childhood Flourishing Index (CFI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Calculate the Childhood Flourishing Index (CFI)
    CFI = rowSums(select(., k6q71_r, k7q85_r, k7q84_r) == 1, na.rm = TRUE),
    
    # Recode CFI into a dichotomous variable
    CFI_dich = case_when(
      CFI == 3 ~ "Flourishing",
      CFI %in% 0:2 ~ "Not Flourishing",
      TRUE ~ NA_character_
    ),
    CFI_dich = factor(CFI_dich, levels = c("Not Flourishing", "Flourishing"))
  )

# Calculate and recode Family Resilience and Connection Index (FRCI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode FRCI items: 1 point for "All of the time" (response = 1), handle missing values
    talkabout_score = ifelse(!is.na(talkabout) & talkabout == 1, 1, 0),
    wktosolve_score = ifelse(!is.na(wktosolve) & wktosolve == 1, 1, 0),
    strengths_score = ifelse(!is.na(strengths) & strengths == 1, 1, 0),
    hopeful_score = ifelse(!is.na(hopeful) & hopeful == 1, 1, 0),
    
    # Recode additional FRCI items: 1 point for "Very well" (response = 1), handle missing values
    k8q21_score = ifelse(!is.na(k8q21) & k8q21 == 1, 1, 0),
    k8q30_score = ifelse(!is.na(k8q30) & k8q30 == 1, 1, 0),
    
    # Calculate total FRCI score by summing all items, ignoring missing values
    FRCI = rowSums(across(c(talkabout_score, wktosolve_score, strengths_score, hopeful_score, k8q21_score, k8q30_score)), na.rm = TRUE),
    
    # Create categorical version of FRCI
    FRCI_cat = factor(
      case_when(
        FRCI <= 1 ~ "0-1",
        FRCI <= 3 ~ "2-3",
        FRCI <= 6 ~ "4-6",
        TRUE ~ NA_character_
      ),
      levels = c("0-1", "2-3", "4-6")
    )
  )
1
ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.
2
ace1 = Hard to Cover Basics Like Food or Housing; ace3 = Child Experienced - Parent or Guardian Divorced; ace4 = Child Experienced - Parent or Guardian Died; ace5 = Child Experienced - Parent or Guardian Time in Jail; ace6 = Child Experienced - Adults Slap, Hit, Kick, Punch Others; ace7 = Child Experienced - Victim of Violence; ace8 = Child Experienced - Lived with Mentally Ill; ace9 = Child Experienced - Lived with Person with Alcohol/Drug Problem; ace10 = Child Experienced - Treated Unfairly Because of Race; ace11 = Child Experienced - Treated Unfairly Because of Health Condition; ace12 = Child Experienced - Treated Unfairly Because of their Sexual Orientation or Gender Identity
3
k2q33b = Anxiety Currently; k2q32b = Depression Currently; k2q34b = Behavior Problems Currently
4
k6q71_r = Show Interest and Curiosity; k7q85_r = Stays Calm and In Control When Challenged; k7q84_r = Works to Finish Tasks Started
5
talkabout = Facing Problems - How Often Talk Together; wktosolve = Facing Problems - How Often Work Together; strengths = Facing Problems - How Often Draw on Strengths; k8q21 = Share Ideas or Talk About Things That Matter; k8q30 = How Well Handling Demands of Raising Children

Frequencies & Descriptive Statistics


Sample Demographic Distributions

Calculate and visualize the distribution of child age, educational attainment, sex, socioeconomic status, and race among neurodivergent participants. Create a raincloud plot for age distribution alongside a summary statistics table, and generate bar plots for key categorical variables, ensuring all visuals are clear and informative.

Expand Code
# Calculate descriptives of age for neurodivergent participants
neurodiv_age_summary <- neurodivergent_data %>%
  summarise(
    Mean = mean(sc_age_years, na.rm = TRUE),
    n = n(),
    SE = sd(sc_age_years, na.rm = TRUE) / sqrt(n()),
    SD = sd(sc_age_years, na.rm = TRUE),
    Minimum = min(sc_age_years, na.rm = TRUE),
    Q1 = quantile(sc_age_years, 0.25, na.rm = TRUE),
    Median = median(sc_age_years, na.rm = TRUE),
    Q3 = quantile(sc_age_years, 0.75, na.rm = TRUE),
    Maximum = max(sc_age_years, na.rm = TRUE),
    IQR = quantile(sc_age_years, 0.75, na.rm = TRUE) - quantile(sc_age_years, 0.25, na.rm = TRUE),
    .groups = 'drop'
  )

# Create raincloud plot displaying child age descriptives
raincloud_plot <- neurodivergent_data %>%
  ggplot(aes(x = 1, y = sc_age_years)) +
  PupillometryR::geom_flat_violin(
    adjust = 1.5, trim = FALSE, alpha = 0.5, color = NA, fill = "deeppink", 
    position = position_nudge(x = 0.2)
  ) +
  geom_boxplot(
    width = 0.12, outlier.color = "black", colour = "black", fill = NA, alpha = 0.5, 
    position = position_nudge(x = 0.2)
  ) +
  stat_dots(
    dotsize = 0.1, side = "left", justification = 1, binwidth = 1, alpha = 0.5, 
    colour = "deeppink", fill = "deeppink"
  ) +
  geom_point(data = neurodiv_age_summary, aes(x = 1, y = Mean), shape = 18, 
             color = "deeppink", size = 3, position = position_nudge(x = 0.2)) +
  geom_errorbar(data = neurodiv_age_summary, aes(x = 1, y = Mean, ymin = Mean - SE, ymax = Mean + SE), 
                width = 0.1, color = "deeppink", position = position_nudge(x = 0.2)) +
  labs(title = "Age Distribution of Neurodivergent Participants", y = "Age (Years)", x = "") +
  scale_y_continuous(breaks = seq(5, 18, by = 1), limits = c(5, 18)) +
  coord_flip() +
  theme_minimal() +
  theme(axis.text.y = element_blank(), legend.position = "none")

# Create summary statistics table as a data frame and transpose it
summary_df <- neurodiv_age_summary %>%
  select(Mean, SD, SE, Median, Minimum, Maximum, IQR, Q1, Q3) %>%
  t() %>%
  as.data.frame() %>%
  mutate_all(~ round(., 2))

# Update row names to meaningful labels
colnames(summary_df) <- "Value"
summary_df$Statistic <- rownames(summary_df)
summary_df <- summary_df %>%
  select(Statistic, Value)

# Convert the transposed summary statistics data frame to a table using gridExtra
summary_table <- gridExtra::tableGrob(summary_df, rows = NULL)

# Combine the plot and summary table using gridExtra
gridExtra::grid.arrange(raincloud_plot, summary_table, nrow = 1)

Expand Code
# Bar plot for educational attainment counts categories
ggplot(neurodivergent_data %>% filter(!is.na(caregiver_education)), aes(x = caregiver_education)) +
  geom_bar(aes(fill = caregiver_education), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  labs(
    title = "Distribution of Adults' Highest Educational Attainment",
    x = "Educational Attainment",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code
# Bar plot for sex counts distribution
ggplot(neurodivergent_data, aes(x = sc_sex_cat)) +
  geom_bar(aes(fill = sc_sex_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  scale_fill_manual(values = c("Male" = "lightblue", "Female" = "pink")) +
  labs(
    title = "Sex Distribution of Neurodivergent Participants",
    x = "Sex",
    y = "Count"
  ) +
  theme_minimal()

Expand Code
# Bar plot for SES counts categories
ggplot(neurodivergent_data, aes(x = fpl_category)) +
  geom_bar(aes(fill = fpl_category), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    title = "Distribution of Socioeconomic Status (FPL)",
    x = "Federal Poverty Level Groups",
    y = "Count"
  ) +
  theme_minimal()

Expand Code
# Bar plot for race counts categories
ggplot(neurodivergent_data, aes(x = sc_race_cat)) +
  geom_bar(aes(fill = sc_race_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.3, size = 3.5) +
  labs(
    title = "Race Distribution of Neurodivergent Participants",
    x = "Race",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code
# Remove intermediate data frames from the global environment
rm(neurodiv_age_summary, raincloud_plot, summary_df, summary_table)

Distribution & Cross-Tabulations of ACEs

Create a bar plot to visualize the distribution of children by ACE counts. Conduct cross-tabulations between ACE counts and psychopathology/flourishing variables, as well as between psychopathology outcomes and the flourishing variable, displaying the results in table format.

Expand Code
# Bar plot of ACE counts categories
ggplot(neurodivergent_data, aes(x = factor(ace_counts), fill = factor(ace_counts))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "ACE Counts", y = "Number of Children", fill = "ACE Counts",
       title = "Distribution of Children by ACE Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))
1
Cross-tabulations (or contingency tables) summarize the relationship between two or more categorical variables. In this case, I’m exploring how the mental health conditions (Anxiety, Depression, Behavioral Issues) and Child Flourishing Index (CFI) are distributed across different levels of ACEs.

Expand Code
# Cross-tabulation between psychopathology/flourishing variables and ACEs count variable
crosstable(neurodivergent_data, c(Anxiety, Depression, Behavioral_Problems, CFI_dich), by = ace_counts, showNA = "no", percent_digits = 1, total = "both") %>%
  as_flextable()

# Cross-tabulation between psychopathology outcomes and flourishing variable
crosstable(neurodivergent_data, c(Anxiety, Depression, Behavioral_Problems), by = CFI_dich, showNA = "no", percent_digits = 1, total = "both") %>%
  as_flextable()

label

variable

ace_counts

Total

0 ACEs

1 ACE

2-3 ACEs

4+ ACEs

Anxiety

No

383 (27.0%)

375 (26.5%)

407 (28.7%)

252 (17.8%)

1417 (7.7%)

Yes

3389 (20.0%)

4366 (25.8%)

5180 (30.6%)

3992 (23.6%)

16927 (92.3%)

Total

3772 (20.6%)

4741 (25.8%)

5587 (30.5%)

4244 (23.1%)

18344 (100.0%)

Depression

No

310 (20.3%)

349 (22.8%)

486 (31.8%)

383 (25.1%)

1528 (16.6%)

Yes

1004 (13.1%)

1475 (19.3%)

2553 (33.4%)

2619 (34.2%)

7651 (83.4%)

Total

1314 (14.3%)

1824 (19.9%)

3039 (33.1%)

3002 (32.7%)

9179 (100.0%)

Behavioral_Problems

No

592 (21.7%)

692 (25.4%)

875 (32.1%)

569 (20.9%)

2728 (14.5%)

Yes

3052 (18.9%)

4234 (26.3%)

4939 (30.6%)

3890 (24.1%)

16115 (85.5%)

Total

3644 (19.3%)

4926 (26.1%)

5814 (30.9%)

4459 (23.7%)

18843 (100.0%)

CFI_dich

Not Flourishing

13914 (27.4%)

14990 (29.6%)

13868 (27.3%)

7936 (15.7%)

50708 (93.0%)

Flourishing

1530 (40.1%)

1239 (32.4%)

808 (21.2%)

243 (6.4%)

3820 (7.0%)

Total

15444 (28.3%)

16229 (29.8%)

14676 (26.9%)

8179 (15.0%)

54528 (100.0%)

label

variable

CFI_dich

Total

Not Flourishing

Flourishing

Anxiety

No

1332 (94.0%)

85 (6.0%)

1417 (7.7%)

Yes

16579 (97.9%)

348 (2.1%)

16927 (92.3%)

Total

17911 (97.6%)

433 (2.4%)

18344 (100.0%)

Depression

No

1477 (96.7%)

51 (3.3%)

1528 (16.6%)

Yes

7530 (98.4%)

121 (1.6%)

7651 (83.4%)

Total

9007 (98.1%)

172 (1.9%)

9179 (100.0%)

Behavioral_Problems

No

2633 (96.5%)

95 (3.5%)

2728 (14.5%)

Yes

15943 (98.9%)

172 (1.1%)

16115 (85.5%)

Total

18576 (98.6%)

267 (1.4%)

18843 (100.0%)

Inferential Statistics


Interaction Models

Fit three logistic regression models with survey weights to examine the moderating effect of Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and psychopathology outcomes (anxiety, depression, behavioral problems). Extract residuals, check heteroscedasticity using Breusch-Pagan tests, calculate Cook’s distance for influential observations, and obtain robust standard errors. Create summary tables with robust standard errors and False Discovery Rate (FDR) adjusted p-values for each model.

  1. The anxiety model did not yield significant interactions, so we do not proceed with stratified post-hoc analyses. Instead, we will run main effect models without interaction terms.
  2. The depression model did not yield significant interactions, so we do not proceed with stratified post-hoc analyses. Instead, we will run main effect models without interaction terms.
  3. The behavioral problems model did not yield significant interactions, so we do not proceed stratified with post-hoc analyses. Instead, we will run main effect models without interaction terms.

Main Effect Models

Fit three logistic regression models with survey weights to examine the main effects of ACE counts and the Child Flourishing Index (CFI_dich) on psychopathology outcomes (anxiety, depression, behavioral problems). Extract residuals, check heteroscedasticity using Breusch-Pagan tests, calculate Cook’s distance for influential observations, and obtain robust standard errors. Create summary tables with robust standard errors and False Discovery Rate (FDR) adjusted p-values for each model.

Session Information


To enhance reproducibility, details about the working environment used for these analyses can be found below.

Expand Code
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS 15.1.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggdist_3.3.2     lubridate_1.9.3  forcats_1.0.0    stringr_1.5.1   
 [5] purrr_1.0.2      tibble_3.2.1     tidyverse_2.0.0  ggpubr_0.6.0    
 [9] crosstable_0.8.1 psych_2.4.6.26   gt_0.11.1        ggeffects_1.7.0 
[13] gtsummary_2.0.3  sandwich_3.1-0   lmtest_0.9-40    zoo_1.8-12      
[17] survey_4.4-2     survival_3.6-4   Matrix_1.7-1     sjmisc_2.8.10   
[21] sjlabelled_1.2.0 ggplot2_3.5.1    kableExtra_1.4.0 knitr_1.48      
[25] readr_2.1.5      tidyr_1.3.1      dplyr_1.1.4      pacman_0.5.1    

loaded via a namespace (and not attached):
 [1] DBI_1.2.3               mnormt_2.1.1            gridExtra_2.3          
 [4] rlang_1.1.4             magrittr_2.0.3          compiler_4.4.1         
 [7] systemfonts_1.1.0       vctrs_0.6.5             quadprog_1.5-8         
[10] pkgconfig_2.0.3         crayon_1.5.3            fastmap_1.2.0          
[13] backports_1.5.0         labeling_0.4.3          utf8_1.2.4             
[16] PupillometryR_0.0.5     rmarkdown_2.27          markdown_1.13          
[19] tzdb_0.4.0              haven_2.5.4             ragg_1.3.2             
[22] bit_4.0.5               xfun_0.45               labelled_2.13.0        
[25] jsonlite_1.8.8          uuid_1.2-0              broom_1.0.6            
[28] parallel_4.4.1          R6_2.5.1                stringi_1.8.4          
[31] car_3.1-2               Rcpp_1.0.13-1           base64enc_0.1-3        
[34] splines_4.4.1           timechange_0.3.0        tidyselect_1.2.1       
[37] rstudioapi_0.16.0       abind_1.4-5             yaml_2.3.9             
[40] lattice_0.22-6          withr_3.0.2             flextable_0.9.7        
[43] askpass_1.2.0           evaluate_0.24.0         zip_2.3.1              
[46] xml2_1.3.6              pillar_1.9.0            carData_3.0-5          
[49] checkmate_2.3.1         insight_0.20.5          distributional_0.5.0   
[52] generics_0.1.3          vroom_1.6.5             hms_1.1.3              
[55] commonmark_1.9.1        munsell_0.5.1           scales_1.3.0           
[58] glue_1.8.0              gdtools_0.4.0           tools_4.4.1            
[61] data.table_1.15.4       ggsignif_0.6.4          mitools_2.4            
[64] cards_0.3.0             colorspace_2.1-0        nlme_3.1-164           
[67] cli_3.6.3               textshaping_0.4.0       officer_0.6.7          
[70] fontBitstreamVera_0.1.1 fansi_1.0.6             broom.helpers_1.17.0   
[73] viridisLite_0.4.2       svglite_2.1.3           gtable_0.3.5           
[76] rstatix_0.7.2           sass_0.4.9              digest_0.6.36          
[79] fontquiver_0.2.1        htmlwidgets_1.6.4       farver_2.1.2           
[82] htmltools_0.5.8.1       lifecycle_1.0.4         fontLiberation_0.1.0   
[85] openssl_2.2.0           bit64_4.0.5